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NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


TECHNICAL REPORT R-178 

THEORETICAL DETERMINATION OF THE BOUNDARY 
AND DISTORTION OF THE GEOMAGNETIC FIELD 
IN A STEADY SOLAR WIND 
By Benjamin R. Briggs and John R. Spreiter 

SUMMARY 


An approximate formulation of the steady-state Chapman -Ferraro problem, 
given recently by Davis and Beard, is used to calculate the coordinates of the 
complete boundary of the geomagnetic field. Field lines are then computed in the 
magnetic meridian plane containing the free -stream direction of the solar wind, 
taking into account the distorting effects of currents flowing in the boundary. 
Numerical results are given for the case where the geomagnetic dipole axis is 
perpendicular to the direction of the solar wind. 

INTRODUCTION 


This paper is concerned with the theoretical determination of the steady - 
state shape and location of the interface, or boundary, between the magnetosphere 
and the solar wind, and of field lines in the confined geomagnetic field. Recent 
accounts of the basic physics, and of the formulation of unsteady as well as 
steady-state mathematical models relating to this problem, will be found in 
references 1 through 5 • 

The present work is an extension of that reported in references 6 through 13 
wherein the three-dimensional Chapman -Ferraro problem is simplified by the intro- 
duction of an approximate boundary condition. The coordinates of the entire 
boundary of the magnetosphere are calculated for the case in which the dipole 
axis is perpendicular to the free -stream direction. A number of field lines are 
then computed in the magnetic meridian plane containing the free -stream direction, 
taking into account the effect of electric currents in the boundary. 

COMPUTATION OF THE COORDINATES OF THE BOUNDARY 
FORMULATION OF THE PROBLEM 

The determination of the shape of the boundary of the geomagnetic field, and 
the total magnetic field B inside it, involves the solution of the magnetic 
field equations div B = 0 and curl B = 0. The earth's magnetic field is repre- 
sented by a three-dimensional dipole singularity at the origin (the center of the 



earth). The normal component of B vanishes at the Boundary, and Dungey (ref. 3 ) 
has shorn that the total (tangential) field at the Boundary, B g , may Be expressed 
mathematically By the relation 


^ o o 

= 2 mnv^ cos y 


(1) 


The quantities m, n, and v are mass, number density, and velocity of the posi- 
tive ions in the solar stream, and f is the angle Between the free -stream 
velocity vector and an outward normal to the surface. The condition cos ^ < 0 
must hold on the Boundary. 

It is a property of the Boundary-value problem described above that the 
field B can vanish only at isolated points on the Boundary. These points are 
designated neutral points. It follows from equation (l) that cos ^ vanishes, 
and the Boundary is therefore parallel to the stream, at these points. 

Beard (ref. 6 ) dropped the condition that the normal component of B 
vanishes at the Boundary and replaced it with the approximate condition that 
B s = 2B-fc, where B-^ is the tangential component of the geomagnetic dipole field 
Bp at the Boundary. The closely related approximation, suggested By Ferraro 
(ref. 9 ), 


B s = 2 fBt ( 2 ) 

where f is a constant, was used by Spreiter and Briggs (refs. 10, 11, and 12) 
and is employed in the present work. 


The differential equation that defines the shape of the Boundary of the 
magnetosphere according to this approximation is obtained By substitution of 
equation (2) into equation (l). It is, following Davis and Beard (ref. 13 ), 
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The variables r , 0, and 9 are spherical coordinates that are fixed with respect 
to the geomagnetic dipole axis (see fig. 1 ). The quantity Mp, the magnetic 
moment of the dipole , is equal to a 3 Bp Q , where a represents the radius of the 
earth, and Bp Q is the magnitude of Bp at the magnetic equator on the surface 
of the earth. The quantities in equation ( 4 ) are in cgs units, and numerical 
values for ay Bp Q , and m are 6.37x10 s cm, 0-312 gauss, and 1 . 67x10 gram 
(for protons), respectively. 

The quantity r G is the geocentric distance along the sun-earth line to the 
boundary of the geomagnetic field. Representative quiet time values for v and n 
are 500 km/sec and 2.5 protons/cm 3 , according to preliminary analysis of data 
from Mariner II (Neugebauer and Snyder, ref. 14 ). These lead to a value for r Q 
of 9*6 earth radii for f = 1 . Axford (ref. 15 ) (see also Kellogg, ref. l6) has 
suggested that 'if the weak interplanetary magnetic field is taken into account, a 
collision -free shock wave may exist in the solar stream on the sunward side of 
the magnetosphere. A consequence of the presence of such a shock wave is that 
the dimensions of the boundary are greater by a factor of 2 1/6 than those indi- 
cated for the present model. Thus the value for r 0 would be about 10.8 earth 
radii in the example given above. 

If attention is confined to the plane cp = ±tc/2(x = 0) where the derivative 
dp/Scp is zero by consideration of symmetry, equation ( 3 ) reduces to an ordinary 
differential equation that can be solved analytically . The trace of the boundary 
in this plane is illustrated in figure 2(a). The front portion is circular, 
p = 1, and the upper portion is defined by the relation 

p eos e = (3/2 2/3 )p 3 /(l + p 3 ) 


The point labeled N is the point of intersection of the two segments of this 
trace. It corresponds to a neutral point of the solution of the exact Chapman - 
Ferraro problem, since the magnetic field indicated by equation (2) is directed 
in opposite directions on either side of these points. All field lines in the 
boundary converge to the neutral points, turn sharply, and extend to the earth. 
The magnitude of the magnetic field should be zero at the neutral points in the 
boundary, although this condition is not fully attained in the approximate 
solution given here. 

The trace in the plane 0 = jt/2(z = 0), the magnetic equatorial plane, is 
obtained similarly by solving numerically the ordinary differential equation to 
which equation ( 3 ) reduces in this case. Tabulations of this result are given by 
Beard (ref. 6) and by Spreiter and Briggs (ref. 10). It is illustrated in 
figure 2(b). 


Notice that the 
equation ( 3 ) • Thus , 


derivatives 8p/c )Q and Sp/8cp appear to the second power in 
it may be rewritten conveniently in either the form 
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The quadratic equations, equations (5) and (6), may be solved algebraically 
for the indicated variables, thus leading to the differential equations 
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It should be noted that equations (9) and (10) are each, in actuality, two 
differential equations, due to the appearance of both the plus and minus signs 
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preceding the radical terms. The trace in the plane cp = ±jt/2, discussed previ- 
ously (see fig. 2(a)), is given hy solutions of the two ordinary differential 
equations to which equation (9) reduces in this case. The circular portion is a 
solution to the equation with the upper sign, and the remaining portion of this 
trace is a solution to the equation with the lower sign. As noted previously, 
the two portions intersect at a point forward of the polar axis that corresponds 
to a neutral point. The magnetic field in this trace of the boundary is pro- 
portional to [sin 0 + 2(cos 0/p) (Bp/50 )] /p 3 , and it is evident that the alter- 
native signs are related to the reversal of the direction of the field vector at 
the point of intersection of the two segments of this trace. 

It is to be expected, then, that the complete approximate boundary will be 
composed of two surfaces that intersect, or in some manner merge, in the vicinity 
of the point over the pole, and that one of these surfaces will be a solution of 
equation (9) (or (10)) with the upper sign, and that the other surface will be a 
solution with the lower sign. The curve on the boundary along which the signs 
must be switched cannot be determined simply, so both surfaces must be computed 
independently. The boundary is then considered to be the exterior part of the 
two intersecting surfaces. 


NUMERICAL INTEGRATION OF THE DIFFERENTIAL EQUATIONS 


The solutions to equation^ (9) and (10) that represent the magnetosphere 
boundary are determined by an integration technique based on Euler ! s method for 
solving ordinary differential equations. (See, e.g., Kunz, ref. 17*) The appli- 
cation of this technique to the present problem will be described specifically 
for equation (9)* It is to be understood, however, that the description applies 
as well for the solution of equation (10) except that the roles of the variables 
0 and cp are interchanged. 

Equation (9) is of the general form 

dp/d0 = F[p, 9, cp, (dp/dcp)] (ll) 

The integration of this equation commences at a known trace of the desired solu- 
tion in some surface 0 = constant, for example, the plane 0 = rc/2 in the 
present problem. Values for p are given at increments Ap along this trace, 
and the derivative bp/btp is computed by numerical differentiation. The deriv- 
ative bp/bd is now readily determined by means of equation (ll). A step A0, 
from the surface 0=0^ to the surface 0 = 0j_+i, is taken by inserting the 
values for p^, 0^, cp^, and (bp/bQ)^ into the linear extrapolation formula 

p i+i = Pj_ + A0(dp/d0 )^ (12) 

The derivatives (dp/dcp)i+i and (3p/d0)i+i are now calculated by the methods 
stated above, and refined values for Pj_ +1 are obtained by substitution into the 
averaging extrapolation formula 
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p i+i = p i + (^/ 2 )[(Sp/3e) i + (dp/de ) 1+1 ] 


(13) 


New values for the derivatives of P-j_ +1 are computed as before, and the 
extrapolation -integration process is repeated for succeeding steps A 0. 

Attention should be called to the fact that the coefficients Aj. and A 2 in 
equation (9) vanish for 0 - 90 ° . Thus, equation (9) is indeterminate at the 
trace in the plane 0 = 90°, and tends to be sensitive to small errors near this 
trace* The derivative 8p/8cp, equation (10), is zero in the plane (p = ±90°, by 
consideration of symmetry. The numerator of equation (10) should therefore 
vanish at cp = ±90°. Small errors in the vicinity of the trace in the plane 
cp = ±90° may lead to negative values for the discriminant B 2 2 -4 B1B3, however, 
thus limiting the use of equation (9) near these traces. Both equations are 
poorly conditioned near the polar axis, in that small errors may lead to negative 
numbers under square root signs. It is due to these complications that integra- 
tions could not be started in the polar region of the trace in the plane cp = 90°, 
nor from the segment cp > l80° of the trace in the plane 0 = 90 ° . Thus, the 
boundary could not be determined by single integrations commencing on one or the 
other of the planes of symmetry, but was calculated in several pieces, as 
described in the next section of this paper. 

The numerical computations have been carried out with the use of an IBM 7090 
computer in the region 90° < cp < 270°, 0 <0 < 90° • The machine procedure 
employed for the taking of derivatives numerically is SHARE Subroutine CL SMB 3, 
Distribution no. 331, which has been converted for FORTRAN use. This program 
involves a 7~P°i rr t polynomial smoothing process, followed by a three -point dif- 
ferentiation process, and it has been found to give good results in quite general 
application. 


RESULTS AND DISCUSSION 


The sunward, nearly hemispherical, portion of the boundary was calculated 
with the use of equation (9), starting on the segment 90° < cp < l80° of the trace 
in the plane 0 = 90° (see fig* 2(b)). The integration was performed in the 
direction of decreasing 0. The upper sign was chosen to precede the radical 
term in the differential equation, by analogy with this choice in the computation 
of the circular part of the trace in the plane cp = 90°. The increments A0 and 
Ap were 2.5° and 5°j> respectively. At 0 = 20° the integration was halted 
because of the occurrence of negative values of the discriminant A 2 2 -kA x A 3 for 
cp in the range 90° < cp < 135°* Two additional integrations were carried out 
starting from the segments 135° < cp < l80° and 150° < cp < l8o° of the trace in 
the plane 0 = 90°. The increments A0 and Ap were taken to be 2° and 2.5°; 
respectively, in these two computations. The first of these was stopped at 
0 = l8°, and the second at 0 = 14°, because of the occurrence, again, of 
negative values for the above indicated expressions. 

It was noted earlier that the trace in the plane cp = 90° consists of two 
intersecting curves obtained with the use of both of the differential equations 
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Implicit in equation (10). The integrations discussed in the foregoing para- 
graph might possibly have been continued to the point over the pole if it were 
known at precisely what point in each plane cp = constant (e.g., the point N in 
the plane cp = 90 °, shown in figure 2(a)) to switch to the other of the two dif- 
ferential equations . The curve in the surface along which the switching takes 
place cannot be determined by any simple technique,, however, so the polar portion 
of the boundary must be calculated by other means. There remains, also, the com- 
putation of the coordinates in the region l8o° < cp < 270 °. 

It was found that the trace previously obtained in the range 14° < Q < 90° 
in the plane cp = l 80 ° can be joined smoothly to the point p = 2 1/3 on the polar 
axis by means of a curve defined by the simple relation 

p = 2 1/3 + Ke 2 (14) 

The constant K is evaluated by substitution into this relation of the value for 
p at some point, say 0 = 20 °, on the previously determined portion of the trace 
in this plane. Equation (10) was used to calculate the coordinates of the bound- 
ary in the region l 80 ° < cp < 27 0 ° commencing with the trace in plane cp = l 80 ° as 
discussed above. The negative sign is chosen to precede the radical term, by 
analogy with the simpler computation of the trace in the plane 0 = 90°. (See, 
e.g., ref. 6 , 10 , or 11 .) 

The integration was performed in the range 0 < G < 90° in the direction of 
increasing cp, with the increments A0 = 2° and Zkp = 5°* The process yielded 
results up to the plane cp = 255° • A second computation was performed with £6 
taken to be 5°- These results agree to at least three significant figures with 
the first integration. A third integration was carried out in the range 
0 < 0 < 20°, where A0 and Ap were taken to be 2° and 5°, respectively. The 
integration produced results up to the plane cp = 265 ° which agree to at least 
three significant figures with the overlapping computations previously discussed. 

The coordinates of the boundary in the region 90° < cp < l80°, near the polar 
axis, were computed by integrations from an assumed trace in the plane cp = 135 °* 
The starting trace was chosen by trial and error such that the calculated traces 
near the plane cp = 90 ° approached the known trace in the plane cp = 90 ° as 
closely and smoothly as possible. Equation (10) was used here, with the upper 
sign. The starting trace was given in the range 0 < 6 < 20°, and the increments 
A G and Acp were 2° and 1°, respectively. The computations were stopped at the 
plane cp = 99 0 in the approach to the plane cp = 90 ° j and at the plane cp = l 68 ° 
in the approach to the plane cp = l80°. The resulting surface merges smoothly 
into the lower nearly hemispherical portion, as previously calculated, for cp 
greater than about 150°. The values for p at the plane cp = 99° were extrap- 
olated linearly, by means of equation ( 12 ), over an increment Acp = - 9 0 to the 
plane cp = 90° • These results agree to three significant figures with the known 
analytic solution in this plane. 

The results of these integrations are shown in figure 3 in the form of 
traces in the planes cp = 90°, 105°, 120°, . . . , 270°. The traces in the 
planes cp = 90 ° and 27 0 ° are the analytic solutions shown in figure 2 (a), and the 
point labeled N corresponds to the neutral point in the exact solution to the 


7 



Chapman -Ferraro problem. Notice that in the region 90° < cp < 1500 the polar 
portion of the boundary intersects the main sunward portion with discontinuous 
slope This discontinuity unquestionably represents a failure of the pres- 

ent approximation, and would not exist in the exact solution. It is believed, 
however, that this is a local failure, and that the present approximate results 
should be in substantial agreement with solutions to the exact problem over 
nearly all of the boundary. Numerical values for p as a function of G in 
planes of constant cp are presented in table I. 


FIELD LINES IN THE MAGNETOSPHERE 


A computation of field lines in the magnetic meridian plane containing the 
free -stream direction will be described in the following paragraphs. 


According to the Chapman -Ferraro geomagnetic storm theory, the field B 
must exhibit a dipole singularity at the origin and satisfy the relations 
curl B = 0 and div B = 0. At the boundary, which forms as a result of inter- 
action with the solar wind, the normal component of B must vanish and the 
tangential component must satisfy equation (l). If the boundary has been deter- 
mined so as to satisfy these boundary conditions, then the field B inside the 
boundary can be determined by solution of the above boundary-value problem. 
Alternatively, but equivalently, it may be determined by summing the effects of 
the geomagnetic dipole and the currents flowing in the boundary by means of the 
equation 


? = ?p + 


Jxd 



(15) 


where Bp represents the geomagnetic dipole and d is the vector from the point 
at which the field is to be determined to points on the boundary. The symbol J s 
represents the currents in the boundary, and it is related to the tangential 
field B s by the equation 

tig — B s xn/4jt ( l6 ) 

where n is an outwardly directed unit vector normal to the boundary. The 
integration in equation (15) is to be carried out over the entire boundary. The 
results for B by the two methods described above will be identical, and field 
lines may be computed by solution of the differential equations 


dz _^z_ , s 

ds | B | ' ds | B | 1 () 

where By and B z are the y and z components of B, and s is a running 
variable along the field line. 

The exact condition that the normal component of B at the boundary is zero 
has been replaced, in the present boundary computation, by the approximate condi- 
tion given by equation (2). If, now, J s is evaluated by means of equation (l6), 
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where B s is consistent with equation (2) and the approximately determined 
boundary, then field lines may be determined by use of equations (15) and (17)* 
Alternatively, the calculated boundary may be considered to be exact, and the 
field may then be calculated by solution of the previously stated boundary -value 
problem. Since the shape of the boundary has been determined only approximately, 
the two sets of resulting values for the enclosed magnetic field will be differ- 
ent. The first mentioned method, involving the use of equations (15) and (l 6), 
is the simpler of the two to apply, and will be employed here. 

The constant f is arbitrary in the approximate boundary calculation, and 
its value may be chosen so as to improve the compatibility of the boundary and 
the field lines in local regions. In the present case, f was selected to assure 
that none of the field lines cross the boundary on the sunward side. The value 
so chosen is 0.85* 

The results for several field lines are presented in figure 4. A few field 
lines corresponding to the undistorted dipole are also shown for comparison. 

These computations are plotted in terms of the dimensionless variable p. Fig- 
ure ij- is a universal plot, for the given value of f, for all values of n and v. 
The radius of the earth, p 0 , depends upon the quantity r Q (eq. (4)), however. 

The field lines are labeled, for convenience, with the polar angle Q e at which 
they intersect the earth, under the assumption that r Q is 9-0 earth radii. 

This case, which forms the basis for the drawing of the earth in figure 4, 
occurs, for instance, for values of n and v of 2.5 protons /cm 3 and 500 lon/sec, 
assuming that f = O.85. 

A gross characteristic of the computed field is that geomagnetic dipole 
lines of force are compressed on both the daytime and nighttime sides because of 
the magnetic effect of electric currents in the boundary. Another is that field 
lines that intersect the earth within about 7° of the polar axis on the sunward 
side are swept rearward. As stated previously, the constant f was chosen so 
that none of the lines of force cross the boundary on the sunward side. The 
compatibility between the approximate boundary and the field lines is thus seen 
to be good forward of the neutral point N. The pattern of field lines tends to 
recede from the boundary from the neutral point rearward, however. It would 
appear, therefore, that the height z of the boundary indicated by the present 
approximate results is somewhat greater than would be given by an exact solution 
of the Chapman -Ferraro problem. This conclusion is supported, furthermore, by 
the fact that the exact solution would indicate that the boundary is parallel to 
the direction of the undisturbed solar wind at the neutral points and that the 
maximum height of the boundary infinitely far downstream must be twice the height 
at the neutral point. Thus, a lowering of the point N, consistent with the 
implications of the calculated field lines, would be reflected in a general 
decrease of the height of the entire rear portion of the boundary. It is antici- 
pated, however, that except in the immediate vicinity of the boundary, and far 
downstream from the earth, the present field computations should provide a use- 
ful approximation to the results that may be anticipated when the Chapman -Ferraro 
problem is ultimately solved exactly. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., April 26, 1963 
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TABLE I. - VALUES OP p IN PLANES cp = CONSTANT 


e. 







48 

deg 







deg 

90 

105 

120 

135 

150 

165 

180 

195 

210 

225 

240 

255 

270 

0 

1 

.260 

1.260 

1.260 

1.260 

1.260 

1.260 

1.260 

1.260 

1.260 

1.260 

1.260 

1.260 

1.260 

5 

1 

186 

1.186 

1.193 

1.205 

1.220 

1.235 

1.260 

1.280 

1.299 

1.316 

1.330 

1-339 

1.341 

10 

1 

116 

1.121 

1.137 

1.159 

1.190 

1.227 

1.261 

1.300 

1-340 

1-375 

1.404 

1.425 

1.430 

15 

1 . 

051 

1.057 

1.079 

1.112 

1.155 

1.211 

1.264 

I.322 

I.382 

1.438 

1.484 

1.520 

1.528 

20 

1 . 

000 

1.013 

1.044 

1.085 

1.136 

1.197 

1.269 

1.340 

1.425 

1.504 

1.572 

1.624 

1.639 

25 

1 . 

000 

1.009 

1.034 

1.073 

1.126 

1.194 

I.277 

1-367 

1.469 

1-573 

1.667 

1.739 

1.764 

30 

1 . 

000 

1.008 
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Figure 2.- Concluded. 







Figure 3.- Traces of the boundary of the magnetosphere in planes of constant angle cp. 
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